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ABSTRACT 

HR 8799 is a four planet system that also hosts a debris disk. By numerically inte- 
grating both planets and a planetesimal disk, we find interactions between an exterior 
planetesimal disk and the planets can influence the lifetime of the system. We first 
consider resonant planetary configurations that remained stable for at least 7 Myrs 
sans debris disk. An exterior debris disk with only ~ 1% the mass of the outermost 
planet (approximately a Neptune mass) was sufficiently large enough to pull the sys- 
tem out of resonance after 2 to 6 Myrs. Secondly, we consider configurations which 
are unstable in less than a few hundred thousand years. We find that these can be 
stabilized by a debris disk with a mass of more than ~ 10% that of the outermost 
planet. Our two sets of simulations suggest that estimates of the long term stability 
of a planetary system should take into account the role of the debris disk. 



1 INTRODUCTION 



Photometric surveys co nducted with Kepler space telescope 
ijBorucki k, Kochl l201ll 'l in additi on to radial v elocit y sur- 
veys such as those conducted by IWright et all (2009) have 
indicated that multiple planet systems are common. Nu- 
merical integrations can be used to determine if these sys- 
tems are stable and est imate the time to first orbit cross- 
ing event or collis i on (Gozdzicwski & Migaszcwski 
Reidemeister et all 120091; Fabrvckv fc Murrav-Clavl 



2009 



2010 



Marshall et al.ll2010l ; iMarois et alj|2010h . Any configuration 
that has a short lifetime is considered less likely. Conse- 
quently, integrations can be used to place constraints on 
both the orbital elements and masses of the planets. These 
numerical investigations often neglect planetesimals. How- 
ever, extraso lar planetary systems can harbor planet esimal 
debris disks (|Su et alj|200g| : iMoro-Martm et al.ll2010h . 

HR 879 9 is a 1.5 ± 0.3M soi A5V star found 39.4 ± l.Opc 
from Earth (|Marois et al.l [200^ 1. It ha s at least four plan - 
ets which have been directly imaged l|Marois et alJ feoiOV 
Mass estimates for the planets, even taking into account 
HR 8799's youn g age of 60+*p° M yr determined by var- 
ious techniques (|Marois et al.l 20081), or 301 "^ as part of 
the Columba association ( Marois et al.l l201Cf ) , have values 
of 10 ± 3 M j upi ter for HR 8799c, d, and e and 7±i Mj for 
HR 8799b for an assumption of a 60 Myr age, and masses of 
7tl Mj for HR 8799c, d, e and 5Mj for an assumption of 
a 30 Myr age. Projected separations for the planets of HR 
8799e, d, c and b from the st ar are observed to be 14.5, 24, 
38 and 68 AU, respectively (|Marois et alj|201oh . However, 
this measurement lacks a long baseline helpful for constrain- 
ing the planets' positions. Astrometric measurements with a 
much longer baseline taken with Hubble Space Telescope in 
1998 for the outer three planets HR 8799b, c and d confirm 



reported values f or HR 8799b and add new observations for 
HR 8799c and d <|Soummer et alj|201lf ). 

Dynamical studies of HR 8799 indicate that it is 
likely in a 4:2:1 dual mean motion resonance (MMR). 
In other words, the inner two planets are in a 2:1 
mean motion resonance while the outer pair are also in 
a 2:1 mean motion resonance. This architecture is re- 
quired to explain how the system has remained stable 
over its observed age ( Gozdzicwski & Migaszcwski 2009; 
Fabrvckv fc Murrav-Clavl l2010l: iReidemeister et all 120091 : 
Marshall et alj l20ld : IMarois et all |201CT ). Possible orbital 
configurations determined by numerica l simulations are 
summarized bv lMoro-Martm et al.l (|2010t ). 

HR 8799 also has an inner debris disk ranging from 6-15 
AU, an outer debris disk which is thought to ext end from 
90 to 300 AU and a dusty halo out to ~ 1000 AU (|Su et all 
2009). The specific details of this model will be discussed 
further in subsection 1.2. 

Planetesimal debris disks can influence the long term 
stability of a planetary system. Wi thin the context of the 
'Nice' model jTsiganis et all 120051 ). planets migrate due 
to interactions with planetesimals, and instability occurs 
when two planets cro ss a strong mean motion resonance. 
Prii ommes et alj |2008T l considered systems put in resonance 
by a gas disk, a possible scenario explaining HR 8799's cur- 
rent configuration. However, after the gas disk had been 
depleted, they found that planetary interactions with the 
remnant planetesimal disk tended to remove these systems 
from resonances and induce dynamical instability. 

It is in this context that we examine the role of the de- 
bris disk in affecting the stability of the HR 8799 system. 
First, we investigate if it is possible to delay the onset of 
instability for an initially highly unstable orbital configura- 
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tion. Then we consider a configuration with a long lifetime 
and determine whether a debris disk can cause instability. 



1.1 Resonant Structure of HR 8799 

As discussed in the introducti o n, in iti al solutions tested 
by Gozdziewski fc Migaszewski J2009I); iReidemeister et all 
l|2009f ); lFabrvckv fc Murrav-Clavl (|2oTol ) had suggested that 
HR 8799b, c and d were most likely in 4:2:1 dual mean mo- 
tion resonances. Nearly all orbital configurations that re- 
main stable during a reverse integration for the estimated 
age of the system which also agree with the observed or- 
bital elements and nominal masses for the planets minimally 
called for HR 8799c and d to be in a 2:1 MMR. Furthermore, 
the dual mean motion 4:2:1 resonance followed by the inner 
pair of HR 8799c and d in a 2:1 MMR allow for the largest 
possible range of masses that still produce simulations sta- 
ble over the lifetime of the system. Alternate solutions in- 
cluded a 2:1 MMR among the outer pair HR 8799b and 
c a s well as a few finely tuned solutions. In a later analy- 
sis, iMarshall et alJ (|201Cll ) found that when the three planet 
configuration is placed in the 4:2:1 MMR with low planet ec- 
centricities, HR 8799 could survive the observed age of the 
system and potentially longer. Re-reduction of Hubble space 
telesco pe astrometric data of HR 8799 by ISoummer et al.l 
(2011), originally taken in 1998, robustly supports the previ- 
ous 4d:2c:lb dual MMR hypothesis with results which have 
only small departures from the exact integer period ratios 
in previously published data. 

This previous dynamical work p recedes the d i scover y 
of the innermost planet HR 8799e by Marois et all |201tt) . 
However, simulations bv lMarois et all |2010h " have indicated 
that the extra planet only places more restrictions on the 
possible system architectures. In order for those simulations 
to remain stable over the minimum estimated lifetime of the 
HR 8799 system, the planets e, d and c are required to be 
in a 4:2:1 MMR along with masses on the minimum end 
of the estimated values. Note that the 4:2:1 dual mean mo- 
tion resonance is still the most likely configuration for this 
system to be in if long term stability is desired - the differ- 
ence being which planets are in the mean motion resonance, 
their projected masses due to dynamical considerations, and 
the maximum age of the system which remains stable given 
these orbital elements. 



1.2 Distribution and Total Mass of HR 8799's 
Debris Disk 

The total mass and distribution of the debris found in the 
disk affects the dynamics, migration rate, extent of migra- 
tion, and the smoothness of migration. However, the to- 
tal mass of the inner and outer debris disk along with the 
dust halo is not well known for HR 8799. Estimating the 
total mass in debris within a disk is difficult to compute 
accurately for both observational and theoretical reasons. 
Solid (non-gaseous) matter, which emits continuously in the 
jtm to mm wavelengths pro duces nearly all of the radiation 
|Beckwith fc Sargentlll99r3 ). However, while the total cross 
section of the dust particles is many orders of magnitude 
larger than that of larger objects (m to km sized asteroids, 
comets, and planetesimals) , it is these later objects which 



comprise a majority of the total mass of the debris. These 
objects are not bright enough to individually detect in the 
wavelengths that they emit. But estimating the total mass 
of the debris is important to determining the dyn amics. 

The halo and disk dust has been modelled by ISu et all 

(2009). The halo mass is estimated to be 1.9 x 10" 2 M e with 
a radius of up to 1000 AU. Estimating the dust mass of the 
other disk components is more difficult - particularly for 
the cold outer disk because the inner and outer edges are 
not well known. In the case of the inner edge, it was first 
thought to be at 90 AU from temperature estimates while 
the outer radius had been modelled at 300 AU to account 
for all the observational constraints, giving a total dust mass 
of 1.2 x 10~ M®. The inner disk is quite warm at ~ 150A", 
allowing for a very low dust mass estimate of 1.1 x 1O _6 M0. 
A brief summary of the modelled parameters that were used 
including assumed surface densities, inner and outer radii, 
minimum a nd maximum gr ain size, as well as total mass can 
be found in lSu et all (|2009T ) in their Table 3. 

Recent submilli meter observa t ions h ave added some ad- 
ditional constraints. Hug hes et all (1201 ll ) examined HR 8799 
and its debris disk with the Submillimeter Array at 880/im 
- a wavelength which is suitable for examining larger dust 
grains. Low signal-to-noise prevented a full multi-parameter 
modelling of the dust at this wavelength but a combination 
of the SMA data and spectral energy distributions (SEDs) 
seem to rule out a narrow ring of dust and favor a broad 
outer ring starting with an inner edge at ~ 150 AU. 

Scaling up from these modelled dust masses to a to- 
tal debris disk mass can be difficult. Small dust grains 
are typically evacuated in debris disks by a combination 
of Poynting-Robertson drag and radiation blow-out (de- 
pending on grain size) on very short timescales relative to 
the age of an average debris disk. To constantly replenish 
the dust in a system that is m any millions of years old, a 



collisional cascade is required ([Sa fronov & Zviagina 19691: 
Dohnanvil Il969l: IWilliams fc Wetherilll Il994l: iTanaka et all 



19961 : iKenvon fc Luull 19991 ; IKenvon fc Bromlevll200ll . l2002h . 

It is assumed that the differential number distribution 
of the particles in mass is a power law which takes the form 



or simply 



dn(m) = Am a dm 



dn(a) = Ca 2 3a da 



(1) 



(2) 



in radius. The steady-state solution depends on an a index 
of This parameter has been estimated empirically and 
numerically in the previously mentioned literature and, as- 
suming that a is determined only by the mass-dependence 
of the collision rate and that the mod el is self-similar, ca n 
be shown to be that value analytically (|Tanaka et al.lll996l ). 

To estimate the total mass of the disk, we integrate the 
differential number distribution times the mass of these ob- 
jects at a constant density from the minimum to maximum 
sizes in the cascade 



M T = 



M(a)dn(a) 



(3) 



The constant in our number distributi on can be set by using 
the model of the outer debris disk bv lSu et all <|2009h . This 
model includes the minimum and maximum grain sizes that 
were used in order to create a synthetic SED which was 
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matched to the observed data. With the dust mass estimate 
from that model, and a minimum and maximum grain size 
of 10/im and 1000/xm respectively, the constant C can be 
computed. 

To find the total mass of the disk we repeat the previ- 
ous calculation only with a new a ma x that corresponds to 
the largest radius objects in our collisional cascade. The to- 
tal mass that this integrand will yield is entirely set by the 
upper bound because there are many orders of magnitude 
between the minimum and maximum object sizes. Unfortu- 
nately, determining a max is difficult. 

Assuming the cascade operates for the system age 
IQuillen et aiT|2007l ) estimated a max , assuming an alpha pa- 
rameter of 11/6, that used only observable properties of the 
debris disk. As per their equation 16, 
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where h, the disk aspect ratio, and f, the normal disk opacity 
at wavelength A, are our disk observables. Other parameters 
include hh, A, r, Q*o, tage, and f T , which correspond to 
the stellar mass, observation wavelength, the radii at which 
there is a break in the surface brightness profile, specific 
energy, age, and an uncertainty factor. 

However, HR 8799 has no constraints on the scale height 
because the disk is not r esolved. We therefor e adopt the /3 
Pictoris estimate for h by IQuillen et al.l l|2007h . /3 Pictoris is 
also an A star which has a similar age and mass to HR 8799 
as well as an extended debris disk. Estimates for the dust 
mass around /3 P ictoris have been fou nd to be 7 '.8M Moon, or 
about 0.096M® dHolland et alj|l998h . This is simil ar to the 
OTgM ffi suggested bv lSu et all (feoog ) for HR 8799. ISu et all 
(2009) also notes that the amount of excess emission in the 
HR 8799 disk is sim i lar to that of other debris disks around A 
stars (See lSu et al.l (l2006l Vs study of t he evolution of debris 
disks around A stars.) See table 1 by IQuillen et al.l |2007T ) 
and references therein for mass, age, and other parameters 
for f3 Pictoris. 

We note that our biggest uncertainty is therefore in h 
in the above formula, which goes as This makes even 
small errors in estimates of h have a large impact on a top - 
However, the goal is not to determine the exact mass of the 
disk, only to determine what a reasonable range is. 

To find a top we must also compute the opacity at a 
specific wavelength. The opacity can be measured by modi- 
fying our total dust mass integral. The opacity at a specific 
wavelength is a measure of the fractional area covered by 
particles of radius a (i.e. the opacity depends on the num- 
ber of particles per unit area times the cross sectional area). 
In this way we can relate the number of particles of a specific 
radius to the opacity and total disk surface area. Secondly, 
it is possible to equate the differential number distribution 
to N a by 



dN(a) 
dln(a) 



IQuillen et al.ll2007h . Thirdly, we make use of a relation 
which describes how the opacity scales with the radius of 
the object, 



T(a) = T d \ 

\a d 



3-q 



(6) 



where q = 3.5 is equivalent too = j from above and ad and 
Td are the radius and opacity at a specific radius/wavelength 
IQuillen et al.ll2007T ). Substituting the above three relation- 
ships into our mass computation yields an integral that is 
dependent only on the radius of the dust. Using 1000/im 
dust particles as the largest grains we find that the opacity 
at the specific dust particle radius is 



Td 



M d 1 
A pad 



(7) 



N(a) 



(5) 



where Md is the total mass of the dust, A is the total disk 
area, p is the density of the dust grains, and ad is the radius 
of th e largest dust g rain in the model. Using dust masses 
from ISu et alj |2009l ) as well as a disk area computed from 
the inner and outer edges and a maximum grain size from 
the same, along with a density appropriate for dust grains 
p — 1.5 — 3.0— 2-3, we compute a value for the opacity of 

T d ~ 10- 4 . 

With these computed disk observables along with the 
other known parameters for HR 87999 we find an a top value 
of a t0 p ~ lkm. This a t0 p yields a disk mass estimate of 
Mdisk ~ 15OM0 or about one half of a Jupiter mass. Large 
debris disks are often attributed masses in the region of 50 — 
100M e . 

While this estimate would place the disk mass usable in 
our simulations at of order a Jupiter mass or less, there are 
a number of caveats to our estimates. We note that our esti- 
mate for atop is low compared to other disks bv lQuillen et al.l 
( 2007) which place the radius at values anywhere from a few 
to a few hundred times larger. Using an estimate for a top in 
line with those for /? Pictoris and similar disks would yield 
a significantly more massive disk. Also, we recognize that 
the masses of the planets are extreme, with four planets of 
ten Jupiter masses each. It is not inconceivable that the disk 
may be significantly more massive than those observed pre- 
viously. We also note that planetesimals with a radius of 
a t0 p are the largest objects that contribute to the collisional 
cascade. Significant mass could be found in other, more mas- 
sive objects that would have collision times too long to con- 
tribute to dust production. Finally, we note that the major 
phenomenon discussed in the paper occur at masses of one 
Jupiter mass or less in the debris disk. It is not required for 
us to include much of the more massive disks, however, given 
the difficulty of estimating the mass and the general pecu- 
liarity and size of HR 8799, we have included fairly generous 
debris disk masses. 

A second more generic estimate of the largest s i zed ob - 
jects that can be grown was found bv lCuzzi et all |l993h . 
Using a numerical model which uses the Reynold's averaged 
Navier-Stokes equations with both turbulence and full vis- 
cosity, they found that it was possible to create 10-100km 
sized planetesimals. This would place our disk mass estimate 
of order ~ Mj. 

Finally, we can use our own solar system as a reference 
point. The 'Nice' model required approximately 30-50 
to explain the outward migration and eventual configuration 
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of our own system (jTsiganis et alj|2005h . This corresponds 
to a few tenths of a Jupiter mass of debris. 

Total debris disk masses used in our simulations varied, 
but had measurable effects at a Neptune mass, or about 
17 Mg. This puts those simulated disk masses in a similar 
range to that of the 'Nice' model. Larger masses than the 
'Nice' model are justified by way of the discussion above. 
Given the uncertainties in measuring atop, it is very difficult 
to determine if the largest total masses can be used in our 
simulations. 

The total particle number used in our simulations, 1024, 
is also comparable to the 1000-5000 particles used in the 
'Nice' model simulations. Above we noted that the diame- 
ter of objects at the top end of the collisional distribution is 
anywhere between one to several hundred km. An at op value 
similar to ft Pictoris of 180km is about 6.5 times smaller 
than the diameter of Pluto. a top values near 1 km would be 
approximately one thousand times smaller. This suggests 
that our simulated disks which are made up of at least a 
Neptune mass in debris (those that are massive enough to 
have measurable effects on the lifetime of the system) would 
have planetesimals that are larger than those that are pre- 
dicted by the collisional cascade. Therefore, while our plan- 
etesimal masses are close to those used in the 'Nice' model, 
they could be more massive than those in HR 8799's ac- 
tual disk. This may result in more stochastic interactions 
between the planetesimals and planets. However, the plan- 
etesimal masses suggested by the collisional cascade are a 
distribution which could reasonably involve both smaller or 
larger planetesimals then those which we have suggested. 



2 INTEGRATOR 

All simulations were run with the software package QYM- 
SYM^. QYMSYM is a GPU-accelerated hybrid second 
order symplectic integrator which permits close encoun- 
ters simila r to th e Mercury software package developed by 
IChambersI l|l999l ). Like Mercury, QYMSYM uses a second 
order symplectic integrator to advance the positions of all 
pa rticles in a simulatio n in the typical manner desc ribed 
bv lDuncan et al.l (| 19981 ) and iLevison fc Duncan! (|l994h . Ex- 
ploratory work on symplectic maps for N-bod y integrators 
was elucidated by IWisdom fc Holmanl l|l99lf ). Additional 
analytical work on the fo rmulations of the sympl ectic in- 
tegrator can be foun d by ISaha fc Tremainel (|l992T ) and by 
lYoshidal |l990l. I1993T I. Unlike these first symplectic integra- 
tors but similar to Mercury, QYMSYM flags any particle 
from an integration when it is deemed that it has a close 
approach to another massive particle during a time step. 
These are then integrated separately with an adaptive step- 
size conventional integrator. The criterion for closest ap- 
proach is decided based on the Hill radii of the interact- 
ing objects. QYMSYM u ses the 4 th order Hermi te inte- 
gration scheme detailed in lMakino fc Aarsethl l|l992h rather 
than the Bulirsch-Stoe r algorithm used in Mercury. See 
i|Moore fc Quillenl 120111 ) for more details on the QYMSYM 
integrator. 

We note that due to the nature of our CUDA based 
1 See author's Web site for source code. 



code, it is not possible to run less than one block worth of 
particles, even in the case where we wish a massless debris 
disk. Furthermore, due to memory latencies, low particle 
counts are not run efficiently on the GPU. Currently, we do 
not have the ability to disable the 0(N 2 ) kick computation 
or a way to offload low particle counts to the CPU. This 
is something that could be rectified in future revisions of 
the code. Additionally, because of the compilation process, 
our code is restricted to operating on machines which have 
CPU's attached. 

These two caveats have an effect on the number and 
size of the parameter space that we are able to explore. 

The inefficiency of the GPU at low particle count and 
single block restriction has the impact of reducing the length 
for which we can integrate our debris-less test case. We do 
not want to run multiple integrators. Two hybrid- symplectic 
integrators will not get identical answers unless the collision 
integration method is the same. Rather than use another 
hybrid-symplectic integrator (like Mercury) to run debris- 
less test cases, we continue to use our own code. However, 
this does make it difficult to run simulations with as many 
orbital periods as those possible with integrators which can 
be run with very low particle counts. 

The GPU requirement also makes it more difficult to 
test a wide parameter space. A node must be equipped with 
a GPU in order to be able to run our code. The GPU clusters 
that we have access to are much smaller than the CPU clus- 
ters that are available. This limits the number of possible 
trials. 

Last, we note that our integrator is truly 0(N 2 ) - all 
particles feel the effects of all others. This is unlike many 
integrators available in the field which typically do not com- 
pute planetesimal-planetesimal effects. This makes our sim- 
ulations more precise, but does increase the arithmetic in- 
tensity. 

All simulations were run on either dual or quad core 
Intel Core2 Duo architecture CPU's with either GT200 
of GF100 architecture NVIDIA CPU's, both of which are 
CUDA capable double precision architectures. The code can 
easily be optimized to run on any Linux kernel 2.6+ distri- 
bution with appropriate GPU hardware. 

2.1 Accuracy, integration parameters and 
eccentricity corrections 

Timestep sizes were set to 0.016 (out of a possible 2tt or- 
bit) corresponding to 42 days in simulations including HR 
8799e and 93 days in simulations without. We used smooth- 
ing lengths which correspond to radii at least a few times 
larger than the size of the planets assuming a planet and 
planetesimal density of the order of lg/cm 3 and disregard 
any simulation during which the energy conservation (given 
by AE/E) dips below 1.0 x 10~ 4 . Additionally, we are only 
interested in the time to the onset of first instability, and do 
not need to concern ourselves with larger energy errors that 
may occur after a planet-planet interaction or ejection. 

In a few simulations, very high eccentricities were de- 
tected for planetesimals which had been ejected from the 
system during a close approach. Because our integrator 
explicitly conserves angular momentum through the use 
of f and g functi o ns wh en solving Kepler's equations (see 
M oorc fc Quillenl i|201ll ) for more details), these large ec- 
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centricities and their corresponding high velocities produced 
a drift in the location of the center of momentum. This is 
most likely due to the way a hybrid symplectic integrator 
typically updates the particles positions. If a particle has its 
position updated such that it is now in close proximity to 
another particle, it would feel a large force. Typically this 
force is removed by reverse integration of the Hamiltonian 
and the particle is then moved to a more accurate integra- 
tion regime (the Hermite integrator in our case). However, 
if by chance that particle is updated to be very close to 
another, the assumptions made in the perturbation theory 
used to split the Hamiltonian into components, namely that 
Keplerian motion is dominant over interparticle forces, no 
longer applies. This means that the particle will not be cor- 
rectly reverse integrated. This problem is common to many 
symplectic integrators. 

Once a planetesimal has been ejected with a very high 
eccentricity this phenomenon is more likely to repeat. The 
high eccentricity means that its angle of impact will be 
nearly perpendicular to the motion of the other orbiting 
objects exterior to its current position. This newly large ve- 
locity relative to the size of the collision detection criteria 
makes it more likely for a particle to be updated from out- 
side the collision detection criteria to be in close proximity 
to the colliding particle without being previously removed. 

Including a larger collision detection criterion or reduc- 
ing the timestep size of the simulation could resolve this 
problem, but at greatly increased simulation time. Decreas- 
ing timestep size goes directly with total simulation time. 
Due to the high velocities achieved by the particles, it would 
need to be decreased significantly. Increasing the collision 
detection radius, which is some factor of the Hill radius, 
will force more particles to be offloaded to the significantly 
slower Hermite integration routines. Not only is the Her- 
mite integrator 0(N 2 ) (and therefore suffers from non-linear 
increases in simulation time based on increasing particle 
count) , but the number of particles included in the Hermite 
integrator will increase by a factor that goes with the colli- 
sional volume - even small increases in the collision detection 
criteria can lead to significant increases in particle counts. 
Alternatively, increasing the smoothing length can partially 
alleviate this problem but has the negative effect of smooth- 
ing out many of the important planetesimal/planetesimal 
dynamics. 

To correct for this occasional error without greatly in- 
creasing simulation time, we wrote an additional check that 
could remove planetesimals which had eccentricities above 
an arbitrary threshold or those that collided with the star. 
Additionally, particles which are no longer bound and have 
a high enough semi-major axis are removed. Energy error 
ch ecks are common in c omparable literature (see appendix 
of iRavmond et al1l201ll for an example) as is planetesimal 
removal due to ejection. 

This check will have no effect on simulations in which 
planetesimals are not ejected or ejected with realistic veloc- 
ities and eccentricities that are still close enough to interact 
or collide with another particle. Only simulations which have 
particles with extremely large eccentricities will be quanti- 
tatively altered. As we would expect, in simulations where 
these particles are removed, it appears that the dynamics 
are qualitatively identical, only without the aforementioned 
drift in center of momentum. We also note that the center of 



mass is conserved to a high degree of accuracy and we main- 
tain a satisfactory level of energy conservation throughout 
either type of integration. Conservation of center of mass is 
typically on the order of one part in 10 12 or better. We there- 
fore present the results of the corrected and uncorrected 
simulations together - using uncorrected simulations when 
no drift is detected in their outputs and corrected versions 
elsewhere. 

The above timestep choice and numerical energy 
integration error are compara ble to those reported by 
iFabrvckv fc Murrav-Clavl l|201fj| ). 

2.2 Simulation configurations of HR 8799 

We set up two initial configurations to test in our simula- 
tions detailed in each respective section. Generally, given the 
importance of the 4:2:1 MMR in previous work, our simu- 
lations begin with this planetary configuration and have an 
additional debris disk or alternatively use a somewhat mod- 
ified version of this planetary configuration. Minimally, the 
simulations have the innermost two planets in a 2:1 MMR. 

Three planet (b,c,d) plus debris disk simulations 
were based off of the stable c onfigurations discussed by 
jFabrvckv fc Murrav-ClavboioT ). using planet masses of 10, 
10 and 7 Mj up for planets d, c and b. Four planet (b,c,d,e) 
plus debris disk simulat ions are based off t he predicted or- 
bital elements found by |Marois et al.|[201ol ') and use masses 
of 7, 7, 7 and 5 Mj up for planets e, d, c and b. Our debris disk 
is simply a uniform distribution of 1024 equal mass parti- 
cles. Depending on total disk mass desired, the planetesimal 
mass was modified accordingly. 

This later set of initial positions and masses reflects the 
most current observational results bv lMarois et all ll2010h by 
including the fourth planet and the corresponding reduced 
masses required from the dynamical simulations that were 
run. The three planet configuration, while no longer repre- 
sentative of those more recent observations and simulations, 
is both practical and illuminating for several reasons. 

First, we note that the planet which is not included 
in the three-planet configuration is the most recently dis- 
covered innermost object. Unless its presence destabilizes 
the system, we expect its effect on the migration rate of 
the outermost planet to be small in those simulations due 
to the proximity and mass of the debris disk to the out- 
ermost planet. Additionally, we only use the three planet 
configuration for simulations in which the system is arbi- 
trarily placed in an unstable configuration from the begin- 
ning. The four planet configuration has much smaller re- 
gions of stability. While the addition of the fourth planet 
would make it even easier for us to create an initially un- 
stable configuration which shares orbital elements similar to 
those that are observed, it makes it significantly more un- 
likely for a system to move from an unstable region to a 
stable region vi a orbital migration. La st, we recall that the 
recent work by |Soummer et al.ll201ll ) which reduced Hub- 
ble Space Telescope observations suggests best orbital fits 
for HR 8799b, c,d that coincide with the 1:2:4 mean motion 
resonance. This data agrees with previous work that fits 
only the three planets, but does not agree wi th dynamical 
simu la tions which include d all four planets (jMarois et al.l 
|201Q[ ). |m arois et al.1 l20ld found that the innermost three 
planets, HR 8799c, d,e were the most likely planets to share 
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mean motion resonances with b excluded. There is as of yet 
no consensus on the possible orbital elements or dynamical 
structure of HR 8799 other than the most long-lived systems 
having the innermost three planets in a 4:2:1 mean motion 
resonance or at least having the innermost two planets in a 
2:1 mean motion resonance. Therefore we assume that the 
2:1 mean motion is the primary initial requirement for our 
simulated systems. 

Because we are measuring the effects of planetary mi- 
gration on system stability to determine if all unstable con- 
figurations remain so, it is useful to run the three planet 
system - it has larger regions of stability for the planets to 
migrate into. Due to the computational intensity of the simu- 
lations and their corresponding time requirements, we were 
only able to run on the order of hundreds of simulations 
over a few months rather than tens or hundreds of thou- 
sands. Given this limitation, it proved difficult to migrate 
large numbers of four planet configurations into regions of 
stability. This lower simulation number and corresponding 
low number of stabilized configurations would introduce a 
fine tuning problem to our analysis so we do not discuss 
those in greater detail, instead focusing on the three planet 
configurations for the stabilization through migration sce- 
narios. This fine tuning issue could potentially be resolved 
via more stable initial conditions which would require signif- 
icantly less planet migration to move from regions of insta- 
bility to regions of stability. This issue is discussed in further 
detail in our results section. 

In the three planet configurations the inner disk edge 
has a semi-major axis of a m in = 2.5 and an outer disk edge 
of a max = 7.0, corresponding to separations of 60 and 170 
AU. In the four planet configurations the inner disk edge 
has a semi-major axis of a m i n — 6.14 and an outer disk edge 
of a max = 20.69, corresponding to separations of 90 and 
300 AU. This second set matches the estimated outer debris 
disk's observed inner and outer edges. 

The inner disk edge with reduced semi-major axis for 
the three planet configuration was used to encourage rapid 
planet-planetesimal crossings and therefore rapid migration. 
As mentioned previously, we forced our three planet configu- 
ration to become unstable on very short timescales. Because 
of this, we require rapid migration in order for any outcome 
other than planet-planet scattering to be a possibility. The 
inclusion of a debris disk at a position more coincident with 
that which is observed would be possible with either faster 
or more GPUs. This is because configurations that are sta- 
ble on longer timescales - which allow for reduced amounts 
of planet migration to be necessary in order to stabilize an 
unstable configuration - could be used. We could addition- 
ally begin to experiment with the inclusion of the fourth 
planet as well as keep the outermost planet at its observed 
position rather than moving it in to encourage the system 
to become unstable. We found both the removal of the in- 
nermost planet and movement of the inner disk edge were 
helpful in finding these post-migration stable configurations. 

Other simulation parameters include the Hill fac- 
tor for encount e r det ection and the K function (see 
iMoore fc Quillenl ifcoill )). both of which are set to 2.0. The 
K function is an arbitrary function which weights certain 
part of the Hamiltonian to be integrated in order to allow it 
to be broken up into evolutions operators which are other- 
wise not possible. When the distance between two particles 



is large, the value of K goes to zero while when the distance 
between two particles is small, the value of K goes to 1 (or 
vice versa). When K or (1-K) are multiplie d into the respec- 
tive s eparated Hamiltonians described in (|Moore fc Quillenl 
2011), it prevents either from becoming too large and break- 
ing the perturbation theory used to separate them. K is a 
function of Hill radii, but was created by trial and error. 

The distribution of planetesimal semi-major axes is flat 
with probability independent of a within a m i n and a max - 
The initial eccentricity and inclination distributions were 
chosen using Rayleigh distributions with the mean eccentric- 
ity e equivalent to twice the mean value of the inclination i 
and i = 0.01. The initial orbital angles (mean anomalies, lon- 
gitudes of pericenter and longitudes of the ascending node) 
were randomly chosen. 

A group of 150 simulations with three planets were run, 
58 using the eccentricity correction while 92 without. Disk 
mass was varied from 1.0 x 10 -30 to 10 Mj up , although we 
only present results from a disk mass of 1.0 x 10 -3 Mj up and 
larger here. 18 simulations with four planets were run, half 
using the eccentricity-correction while the other half with- 
out. These simulations differed only in the initial conditions 
of the planetesimals. 



3 SIMULATIONS OF INITIALLY UNSTABLE 
PLANETARY ORBITAL ARCHITECTURES 

Does the addition of a relatively massive planetesimal disk 
allow for an increase in stability timescale? To answer this 
question we created an unstable initial configuration for 
both the three and four planet configurations of HR 8799. 
In either case, this was done by starting with actual ob- 
served orbi tal elements f ound by iFabrvckv fc Murrav-Clavl 
l|2010h and|M arois et all (j2010h and reducing the outermost 
planet's semi-major axis in small increments until the sys- 
tem was unstable on the time scale of thousands of orbits or 
less. In the three planet simulations, this gives an initial con- 
figuration with all three pl anets having orbital elements the 
same as those reported bv IFabrvckv fc Murrav-Clavl (|2010h 
except that the outermost planet's semi-major axis was re- 
duced by 14%. For the four planet simulations, the plan- 
ets' orbital configu ration was identical to those found by 
iMarois et alJ |2010) but the outermost planet has a semi- 
major axis reduced by 8%. This arbitrary and largely un- 
stable configuration was chosen to both reduce simulation 
time as well as allow for pronounced effects by migration. 
Due to the previously mentioned limitation in the number 
of total simulations it is possible to run and the extremely 
limited regions of stability available to the four planet con- 
figuration, we opted to run far more three planet simulations 
than four planet systems and present only that data. 

Simulations both with and without eccentricity correc- 
tions are presented simultaneously and see no significant dif- 
ference between them. 

3.1 Results 

In figure [1] we plot the difference between the stability time 
for each simulation and that of a system lacking a debris 
disk. The stability time is the time to first planet-planet 
encounter. The time difference is plotted as a function of 
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Figure 2. In this figure we show the evolution in semi-major 
axis of an example three planet configuration over time, a) Shows 
the evolution in time of the semi-major axis of the three planets 
over the entire 5 Myr simulation, b) Shows the evolution in time 
of the semi-major axis of the three planets over the first Myr. 
In b) we also include error bars that indicate the maximum and 
minimum separation (given by a(l+e) and a(l-e) respectively) 
of the planets given their eccentricity. This second image more 
clearly shows the migration of the outermost planet planet from 
its initial position at 58.6 AU to approximately 65 AU over the 
first few hundred thousand years. In this example the total disk 
mass was 1.6Mj up and the configuration remained stable over 
the entire 5 Myr simulation. 



disk mass. Simulations are run for a maximum of 5 Myrs. If 
no encounters had taken place in that time, they are plotted 
as upper limits in the figure. 

In figure [2] we plot the evolution in time of the semi- 
major axis of the three planets for a sample simulation. In 
this case the total disk mass simulated was 1.6Mj up . 

Figure [1] shows that at low disk masses the difference 
in lifetime is near 0. Only massive disks with masses of ten 
percent or more of the outermost planet can remain stable 
until the end of the simulation. Given enough mass, migra- 
tion via planetesimal scattering can take a tightly packed 
configuration of HR8799 b, c, d and e, or in the case of fig- 
ures [1] and [2l simply b, c and d, and migrate the outermost 
planet into a more stable configuration. The total migration 
for the outermost planet in system configurations which re- 
mained stable over the 5 Myr simulation ranged from 5 to 15 



AU in semi-major axis depending on total disk mass. This 
migration occurred over a few hundred thousand years. 

If stable regions are small compared to unstable re- 
gions, migration is unlikely to put an initially unstable con- 
figuration into a stable region. Only large migrations could 
significantly increase the stability time. This is consistent 
with what we see; rapid migration caused by a massive 
disk is able to significantly affect the lifetime of our simula- 
tions. Migration is expected to reduce planetary eccentric- 
ities and this could increa se stability l|Lissauer &: Steward 
1 19951 : iFernandez felp|ll996h . We searched for signs that the 
outermost planet's eccentricity decreased during migration. 
In figure [2] we see oscillations in eccentricity of the outer- 
most planet were large throughout the simulation due to the 
proximity and high masses of the inner planets. We therefore 
attribute the increased stability time to the wider separation 
rather than any eccentricity damping. 

As a planet migrates outwards it can cross mean motion 
resonances with other planets. For simulations with short 
lifetimes, we searched for evidence that resonant crossing 
caused instability. We examined any relevant 1 st , 2 nd , 3 rd 
and 4" 1 order resonance possible between the outer planet 
and innermost planet as well as the outer planet and the 
middle planet. We saw no large semi-major axis or eccen- 
tricity jumps caused by resonant crossing so they are un- 
likely to be the cause of instability in the simulations with 
shorter lifetimes. Resonances near the outer planet's posi- 
tion are primarily 2 nd or 3 rd order and cause smalle r eccen- 
tricit y jumps due to the weaker dependence on mass (|Quillenl 
12006( 1 

As discussed in our section on simulation parameters, 
our choice of unstable configurations was done in part to 
reduce the amount of simulation time. With an initial plan- 
etary configuration so far from a region of stability, a large 
amount of migration is required in order to move the planet 
to a more stable region of phase space. Because of this, only 
a massive debris disk can migrate the outermost planet suf- 
ficiently far to increase the lifetime. 

Another concern with this set of simulations is that with 
a planetesimal count of 1024 and a total disk mass of at least 
a Jupiter mass, interactions between p lanetesimals and plan- 
ets would cause stochastic mig ration (IZhou et al.ll2002T i. We 
would have expected such large planetesimal masses to cause 
instability rather than increase the lifetime of the system. 
Therefore, we have no reason to suspect that stochastic be- 
havior is a problem. Because the typical disk mass required 
to stabilize the system is unrealistically large there may be 
no need to resolve the disk more accurately. As mentioned 
previously, an approximate upper limit for debris mass is 
thought to be on the order of 10 _3 Mq, or approximately 
one Jupiter mass. For HR 8799, this is effectively the mini- 
mum amount of mass required for any observable effects in 
our simulations. 



4 SIMULATIONS OF INITIALLY STABLE 
PLANETARY ORBITAL ARCHITECTURES 

In this second set of simulations we examine the impact of 
a planetesimal debris disk on the stability time of a con- 
figuration of planets which was stable for about 150k or- 
bits (or around 7 Myr) sans debris disk. In these Simula- 
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Figure 1. This figure shows the difference in time until a planet-planet encounter for a group of HR 8799-like simulations with a 
debris disk from an identical simulation with no disk mass versus the amount of mass in that particular disk. An intentionally unstable 
solution was used in order to q uickly show the effects of a planetesimal debris disk. All simulations use orbital elements taken from 
iFabrvckv fc Murray-C lay (2010) for the inner two planets (which are in a 2d:lc resonance). The outer planet, HR 8799b, has a semi- 
major axis which is reduced by 14% from approximately 68 to 58.6 AU in order to achieve a planet-planet close approach within a couple 
hundred thousand years. At represents the increase or decrease in time before the first planet-planet encounter. Simulations were halted 
after the first planet-planet encounter or after 5 Myr. 



tions we begin with all four planets situated in an orbital 
configuration based off the current estimated positions by 
M arois et alj (|2010h . This four planet configuration is gener- 
ally extremely unstable unless situated in a 4e:2d:lc MMR. 
As before, smoothing lengths are set to sizes a few times the 
radii of the planet to prevent the unrealistic forces possible 
if two particles form a very tight binary or collide within 
what would have been considered an approximate planetary 
radius. Energy error is typically even lower in these simu- 
lations with a AE/E at 1.0 x 1CT 5 or 1.0 x 10~ 6 due to 
reduced planetesimal masses. We stopped simulations after 
their first planet-planet encounter regardless of the even- 
tual fate of the interacting planets. Stability time in the 
histogram is the time to first planet-planet encounter. 

We began by simulating the four planets with an ini- 
tial or bital configuration comparable to that bv lMarois et al.l 
|2010l ) with massless debris and plotted the resonant argu- 
ments. We see constrained resonant angles which indicate 
the presence of a 2:1 MMR between the inner t wo plan- 
ets, shown in figure El While iMarois et all l|2010T i did not 
directly show the resonant argument to illustrate that HR 
8799d an d e were in resonance, ou r plot is similar to Fig- 
ure 10 bv lFabrvckv fc Murrav-Clavl |2010l ) although for HR 
8799e and d rather than d and c. The integrator confirms 
th at this planeta r y con figuration is in resonance as suggested 
bv lMarois et~al] (feoiol ). 



Matching the exact age found by IMarois et al.l (|2010l ) 
is difficult despite the fact that there is no longer a require- 
ment of planetesimals with mass. This is due to the way 
the integrator operates. Specifically, it is not possible to in- 
tegrate less than a certain number of particles, depending 
on compile time and hardware restrictions. For more de- 
tails on the exact natur e of the code, we refer the reader to 
iMoore fc Quillenl (|201ll ). However, while the control simu- 
lation is not as long lived as those found by IMarois et al.l 
l|20ld ). it is still stable on long time scales (10 5 orbits). 



This simulation is then compared to 18 simulations with 
Neptune mass disks. These 18 simulations are identical in 
terms of planetesimal mass, total disk mass, initial outer 
and inner disk edge, and initial planetary orbital configu- 
ration. However, the planetesimals initial orbital configura- 
tions have been randomized 18 different ways. 



The histogram shown in figure [4] shows the distribution 
in lifetimes of the eighteen simulations each with a Neptune 
mass disk. Figure shows the evolution in time of the semi- 
major axis of all four planets of a sample simulation. The 
error bars indicate the minimum and maximum separations. 
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Figure 3. Here we show the resonant angle for 2:1 mean motion 
resonance for HR 8799d and e as function of time, a) For a sim- 
ulation without a debris disk, b) For a simulation with a debris 
disk that has a mass of l/100th the outer planet. The initial con- 
ditions for the planets were stable for 7 Myrs. In b) the extreme 
values of <j> slowly begin to increase over the length of the sim- 
ulation. This effect is not present in a). The resonant island is 
so small that even a disk mass of 1% the outer planet is capable 
of effecting the systems dynamics. The resonant angle plotted is 
given by <j> e = 2X c i — X e — zo e . 

4.1 Results 

In figure |4] we see that in 16 of 18 simulations the stability 
time has decreased by at least 2 Myrs from the initial 'stable' 
configuration which had its first planet interactions after 7 
Myrs. We see a broad distribution of lifetimes. Even though 
a Neptune mass disk is only about 1% the mass of the outer 
most planet it has an effect on the stability time. 

With a Neptune mass disk each simulated planetesimal 
has a mass of about 8 Plutos. Debris disk models estimate 
the top of the collisional cascade to conta in similar (but 
slightl y smaller) size bodies. For example, IWvatt fc DenfJ 
(2002) argued for a distribution of planetesimals in the Fo- 
malhaut system with a range of sizes between 4 and 1000km. 
As discussed previously, this second value roughly coincides 
with Pluto's diameter within a factor of a few of what 
we used in our simulations. Stochastic perturbations are 
expected in real systems; however they will be somewhat 
smaller than those present in our simulation do to our fac- 
tor of 8 difference in mass of planetesimals. 
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Figure 4. Here we show a histogram of time to the onset of 
instability, which we define as planet interaction or ejection. The 
disk mass is about one Neptune in mass. The planets are initially 
in a configuration that is stable for approximately 7 Myrs. We see 
that the distribution of lifetimes is very broad and many of the 
simulations had lifetimes less than this. The reduction in lifetimes 
was as large as ~ 6 Myrs. 




1 ' ' ' ' ' ' 1 

0.5 1 1.5 2 2.5 3 

t [Myr] 

Figure 5. In this figure we show the evolution in time of the 
semi-major axis of an example 4 planet configuration. This par- 
ticular simulation became unstable in a little over 3 Myr. We 
include error bars that indicate the maximum and minimum sep- 
aration (given by a(l+e) and a(l-e) respectively) of the planets 
given their eccentricity. In all four-planet simulations the total 
disk mass was equivalent to Neptune. 

Figure 3 bv lGozdziewski fc Mi aaszcwskil i|2009t ) demon- 
strates why a planetesimal disk with a total mass around 
that of Neptune can have such a pronounced effect on the 
stability time. In this plot, dark regions indicate highly likely 
configurations of eccentricity and semi-major axis for planet 
d, whereas yellow regions indicate strongly chaotic systems. 
Here we see that semi- major axis changes of less than 0.1 AU 
can move planet d from a stable region to a strongly chaotic 
one. Similarly, minor changes in eccentri c ity ca n also have 
a large effect. iFabrvckv fc Murrav-Clavl (|20ld Vs Figure 7 
shows a similar behavior. Note that a change in current sep- 
aration of planet d by 0.05 in the figure (corresponding to 
semi-major axis change of 1.22 AU) can decrease the stabil- 
ity time by as much as four orders of magnitude. A three 



10 



order of magnitude change in stability time is possible with 
a semi-major axis increase of approximately half an AU. In 
either example the regions of stability have sharp edges be- 
tween stable and chaotic solutions. 

We would expect a disk mass of 1% the mass of a planet 
to be able to migrate a planet roughly 1% of the semi-major 
axis if a majority of its angular momentum is transferred to 
the planet via scattering. If this were the case than a Nep- 
tune mass disk would be able to migrate the outer planet 
(which has a mass of 5Adj up at 68 AU) approximately 0.5 
AU. This is more than enough to vary the lifetime by many 
orders of magnitude. Therefore we looked for the number of 
planetesimals which had become orbit crossing by the time 
each simulation had become unstable. In all 16 simulations, 
at least 15% of the disk mass was orbit crossing. This corre- 
sponds to a minimum of 2.5 M®. By starting a simulation 
near the edges of a region of stability, we see that even lower 
mass debris disks can affect the stability time. 

4.2 Pulling HR 8799's planets out of resonance 

We continuously monitor the planets' orbital elements 
throughout all simulations to see if they are in resonance 
with one another. We do this by noting whether or not the 
resonant angle librates around or n. We looked for 1 st , 
2 nd , 3 rd and 4 th order mean motion resonances that could 
exist near each of the planets' semi-major axis at varying 
times. We did not observe any strong two body resonances 
beyond the 2:1 MMR between e and d, a much weaker 2:1 
MMR between d and c, and no clear resonant angles for b. 
We attribute this to the overwhelmingly large gravitational 
perturbations caused by the massive planets as well as the 
fact that most of the nearby resonances for the planets were 
of higher order. This prevented behavior similar to the 'Nice' 
model during which migration causes resonance crossing and 
corresponding eccentricity increases. Finally, we see that in 
all simulations when the planets e and d are pulled out of 
the 2:1 MMR the system rapidly becomes unstable. 

In figure [3^, we see that the resonant angle for the 2:1 
MMR between HR 8799e and d over the 7 Myr simulation 
remains largely unaffected until the end of the simulation. 
However, in figure [3p we note that the extreme values of 
4> begin to increase over time, leading to e and d moving 
out of resonance and the rapid disruption of the system. In 
this particular example, the system's lifetime is decreased 
by about 50%. It is possible the debris disk is responsible 
for pulling HR 8799 out of resonance and so, reducing its 
lifetime. 



5 CONCLUSION 

In this paper we discuss how a planetesimal debris disk can 
effect the stability of the multiple planet system HR 8799. 
Two questions are considered; whether it is possible to desta- 
bilize a stable planetary configuration with a planetesimal 
debris disk and conversely, whether it is possible to stabilize 
an unstable configuration with a planetesimal debris disk. 
We examine both three planet (b,c,d) and four planet con- 
figurations (b,c,d,e). 

In three planet configurations which were unstable on 
short timescales without a debris disk, only massive debris 



disks (10% the mass of the outermost planet) could cause in- 
creases in stability time. In four planet configurations which 
were stable over long timescales without a debris disk, de- 
bris disks of only a Neptune in mass (1% the outer planet) 
can cause large decreases in lifetimes. We attribute this sen- 
sitivity to the small size of HR 8799's resonant region of 
stability. 

While the amount of disk mass required for system sta- 
bilization in our simulations is unrealistic, it is only required 
to be this large due to the initial conditions. In order to run 
a large number of simulations, a highly unstable configu- 
ration that rapidly had planetary encounters was required. 
Given an initial configuration which is just outside a region 
of stability, it would be possible to cause sufficient migration 
to allow an unstable configuration to become stable with a 
more reasonable disk mass. 

Similarly, while the disk properties required for system 
disruption in our simulations are reasonable in both total 
mass and planetesimal size, the mass distribution for the 
planetesimals is not a true mass distribution. In order to 
keep total particle number down we restrict all planetesimals 
to have an identical 8 Pluto sized mass. We do not believe 
that this had a large effect on the distribution of stability 
times of the system when a disk is included because a realis- 
tic mass distribution could encourage even more stochastic 
migration if there are even a few planetesimals of greater 
mass. 

These results suggest that HR 8799 may be destined 
for eventual planet scattering. The very high masses of the 
planets causes regions of stability to be relatively small and 
makes it possible for low mass debris disks to pull the sys- 
tem out of its mean motion resonances and induce instabil- 
ity. We also see that migration is much more likely to make 
the system unstable than migrate the system to a region of 
stability. It may be the case that the debris disk has already 
been responsible for removing the system from a maximally 
stable region and is currently near the boundary of instabil- 
ity. All of these possibilities tend to induce instability rather 
than stability and suggest that HR 8799 may be headed to- 
wards instability, possibly sooner than would otherwise be 
predicted. 

In this study we have used initial conditions consistent 
with observed positions of the planets when attempting to 
determine whether a stable orbital configuration could be 
made unstable by a debris disk. As we have shown here, a 
low mass disk can cause evolution of the planetary config- 
uration. Consequently, in the past the planets could have 
been in a different configuration. By exploring different ini- 
tial conditions future investigations could explore scenarios 
for the past evolution of HR 8799 that would be consistent 
with its current configuration. For example, if the planetes- 
imal disk is causing the system to move out of resonance, 
in the past it may have been in a more stable region of this 
resonance, instead of on its boundary. Increasing simulation 
length and particle count could also increase the resolution 
of our simulations. 

Alternatively, it appears that similar simulations could 
be used to put upper limits on planetesimal debris disk mass 
if it is assumed that the debris disk has a negligible dynam- 
ical effect. In the case of our simulations of HR 8799, plan- 
etesimal disk masses would need to be much smaller than 
one Neptune in mass. 
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The stability timescale for an HR 8799-like system is 
very sensitive to small alterations in planetary orbital ele- 
ments. This sensitivity implies that the effects of planetesi- 
mal debris disks may not be negligible for systems with such 
small regions of stability. However, dynamical simulations of 
HR 8799 sans planetesimal disk were used to determine the 
resonant structure of the system and constrain observations 
while not including all the required dynamics. The assump- 
tion of long term maximum stability may be erroneous when 
debris disks could allow for migration of planets from cur- 
rently observed stable regions to unstable configurations or 
vice versa. In general, it would appear that using dynami- 
cal models which do not include the effects of planetesimals 
to constrain observations may be unwise when a system is 
known to have such limited regions of stability. 

Additionally, the long term stability of systems with 
lower mass disks could still be important. In these systems 
the planets will typically be smaller and the regions of sta- 
bility would be larger. However, reduced planet mass would 
increase the migration rate. 
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